Binding behavior of receptor binding domain of the SARS-CoV-2 virus and ivermectin

The COVID-19 pandemic, caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), sparked an international debate on effective ways to prevent and treat the virus. Specifically, there were many varying opinions on the use of ivermectin (IVM) throughout the world, with minimal research to support either side. IVM is an FDA-approved antiparasitic drug that was discovered in the 1970s and was found to show antiviral activity. The objective of this study is to examine the binding behavior and rates of association and dissociation between SARS-CoV-2 receptor binding domain (RBD), IVM, and their combination using aminopropylsilane (APS) biosensors as surrogates for the hydrophobic interaction between the viral protein and human angiotensin-converting enzyme 2 (ACE2) receptors to determine the potential of IVM as a repurposed drug for SARS-CoV-2 prevention and treatment. The IVM, RBD, and combination binding kinetics were analyzed using biolayer interferometry (BLI) and validated with multiple in silico techniques including protein–ligand docking, molecular dynamics simulation, molecular mechanics-generalized Born surface area (MM-GBSA), and principal component analysis (PCA). Our results suggest that with increasing IVM concentrations the association rate with the hydrophobic biosensor increases with a simultaneous decrease in dissociation. Significant kinetic changes to RBD, when combined with IVM, were found only at a concentration a thousand times the approved dosage with minimal changes found over a 35-min time period. Our study suggests that IVM is not an effective preventative or treatment method at the currently approved dosage.


Preparation of samples containing SARS-CoV-2 RBD
The materials for cell growth and purification of RBD were prepared following the same procedure by Zhang et al. 38 .RBD was expressed in HEK293 cells obtained from Dr. Jason McLellan, Dept. of Molecular Biosciences, The University of Texas, Austin and the Fc and His-tagged recombinant RBD was purified using the Ni-IMAC FF Sepharose column.Using the prepared RBD, which had an original concentration of 1.822 mg/ml and molecular weight of 58 kDa, and Phosphate-buffered saline (PBS), 1:2, 1:5, and 1:10 dilutions of RBD were prepared at 25 °C then stored at 4 °C in Eppendorf tubes.The 1:5 RBD dilution was then heated to 37 °C in a water bath and mixed with various IVM concentrations.

Preparation of samples containing ivermectin
IVM in powder form (molecular weight of 0.8751 kDa) and PBS were used to create a 1:1 dilution of 380.5 μM concentration at 25 °C.Then a serial dilution was performed with PBS to create 1:10, 1:100, 1:1000, and 1:10,000 dilutions, also at 25 °C, and then stored at the same condition as the RBD samples.

Determination of the basic kinetics for SARS-CoV-2 RBD
To determine and study the basic kinetics of the prepared samples of RBD of SARS-CoV-2 and IVM, BLI was performed on individual and mixed samples.The primary instrument used in this study was the Octet® R4 BLI Label-Free Detection system (Sartorius) and APS biosensors (Fortebio).The APS biosensors were hydrated with PBS in a 96-well black plate (Grenier) for at least 30 min before use.The samples were heated in a water bath at 37 °C for at least 20 min before plating 200 uL of each sample into a separate 96-well black plate.The combination samples were prepared by first adding the RBD into the plate and then the IVM before mixing with the pipette tip.The Octet® R4 system is capable of running four tests at a time, each following the same pathways and within a system set at 37 °C with the sample shake plate on.Each test contains a 60 s baseline step with PBS used as the buffer, a 300 s association step with the sample, and then 300 s for dissociation into PBS.It is important that at least one reference test is performed in each experiment by running the entire test with the PBS buffer.The Octet® CFR software was used to calculate and display all binding, association, and dissociation results.The mean and standard deviation of the kinetic constants were shown on the graphs.

Statistical analysis
The one-way Analysis of Variance (ANOVA) was applied using JMP 17 to determine the statistical significance among three or more groups, p-values less than 0.05 were considered statistically significant.A comparison for all pairs using Tukey-Kramer HSD was also performed to determine the means significantly different from each other when applicable.

Protein-ligand docking verification
The crystal structure of the spike protein receptor binding domain of SARS-CoV-2 in complex with ACE2 was obtained from the RCSB protein data bank (rcsb.org,PDB: 6VW1) at 2.68 Å resolution 39,40 .The RBD (Chain E) was extracted and prepared by the Protein Preparation Wizard in Schrödinger suite 41 .The polar hydrogen atoms were added to amino acid residues with valence errors and one missing side chain was fixed by Prime 42 .The IVM structure was obtained from PubChem 43 and prepared by the LigPrep Wizard in Schrödinger suite 44 .Molecular docking was performed in Glide with the docking grid of size 10 × 10 × 10 Å set around the interaction site of RBD and ACE2 complex 45 .The protein-ligand docking poses and interactions were analyzed in Maestro 46 .

Molecular dynamics
Molecular dynamics (MD) simulation was performed on the RBD and IVM conformation with the highest docking score using Desmond 47 .All simulations were performed in three replicates.An orthorhombic simulation box with minimized volume was generated to construct a solvated system using the System Builder in Desmond.Three commonly used explicit water models with three interaction sites-simple point charge (SPC), extended simple point charge (SPC/E), and transferable intermolecular potential 3P (TIP3P)-were initially built, and MD simulation with each model was performed for 10 ns to compare the effects of water models and determine which model to proceed for the longer simulation.The final MD simulation was built under the TIP3P model.Two chloride ions were used to neutralize the system and 0.1 M sodium chloride was added.Prior to the simulation, the system underwent the standard relaxation protocol to reach equilibration.Three protein-ligand complexes were tested-unbound RBD, and RBD with one and five IVM molecules.The full 100 ns MD simulation was performed with the OPLS3e force field in the isothermal-isobaric (NPT) ensemble at a temperature of 300 K and pressure of 1 atm for the recording interval of 100 ps which generated 1000 frames in total 48 .The protein α-carbon (Cα) root-mean-square deviation (RMSD) values were obtained with the Simulation Interactions Diagram Wizard in Desmond.The molecular mechanics-generalized Born surface area (MM-GBSA) method was used to calculate the binding free energy between the RBD protein and ivermectin over the 100 ns simulation.The MM-GBSA calculation was performed using the thermal_mmgbsa.pyscript in Schrödinger Prime which splits the trajectory file generated in Desmond into individual snapshots and computes the average binding energy of the frames 42 .The chosen calculation method used the VSGB 2.0 dissolution model and the OPLS3e force field 48,49 .The initial and final frame from the MD simulation of the RBD and IVM were extracted and analyzed using the Hydrophobic/philic Surfaces Panel.The hydrophobic and hydrophilic surface area of RBD was analyzed with the cut-off particular potential value (isovalue) of − 6 kcal/mol for hydrophilic regions and − 0.5 kcal/mol for hydrophobic regions.The trajectory file was imported into Visual Molecular Dynamics (VMD) software to calculate the solvent accessible surface area (SASA) and the radius of gyration (Rg) of the RBD protein in the different systems to study the protein stability over time 50 .Principal component analysis (PCA) using the pairwise distance method was performed on the RBD protein Cα to study the system motion during the MD simulation.The python package MDTraj was used to align the trajectory snapshots to the initial frame and construct the eigenvectors 51 .

Binding kinetics of ivermectin
The Octet system determined the molar concentrations for IVM at the 1:1, 1:10, 1:100, 1:1000, and 1:10,000 dilutions to be 380.5 μM, 38.05 μM, 3.805 μM, 0.3805 μM, and 0.03805 μM based on the molecular weight of 0.8751 kDa.The basic kinetics analysis was performed at the average body temperature of 37 °C.A higher binding was shown for the highest IVM concentration, in which the 380.5 μM IVM demonstrated dramatic association and dissociation, while the other concentrations displayed binding at a consistently low level (Fig. 1a).The stronger binding capability of IVM with the highest concentration tested resulted in an average binding layer of 0.646 nm, compared to those of the other concentrations ranging from 0.033 and 0.109 nm by the end of the 300 s association step.This binding occurred throughout the entire association step but at a quicker rate toward the beginning.A similar observation of higher binding with increasing IVM concentration was made by Nappi et al. when IVM was tested against purified heat shock proteins 27 (HSP27) using BLI 52 .Previous study of IVM for the inhibition of the epidermal growth factor receptor (EGFR) also reported similar binding trends but at a lower range 53 .HSP27, EGFR, and ACE2 receptors all contain hydrophobic regions of various sizes which could explain the similarities and differences seen in their interactions with IVM [54][55][56] .The average association constants (ka), which measures the rate of sample binding to the hydrophobic biosensor, are shown in Fig. 1b.The average ka of 380.5 μM, 38.05 μM, 3.805 μM, 0.3805 μM, and 0.03805 μM IVM were 5.71 × 10 1 Ms −1 , 5.75 × 10 7 Ms −1 , 5.60 × 10 8 Ms −1 , 5.62 × 10 6 Ms −1 , and 1.33 × 10 8 Ms −1 , respectively.The 380.5 μM concentration was found to associate with the APS biosensors much weaker.Figure 1b demonstrates significantly different ka among tested IVM concentrations with an increase of five to seven magnitudes between 380.5 μM and the other concentrations (p = 0.05).The average dissociation constants (kdis) were also tested and returned similar results.The kdis values demonstrate how easily the IVM separates from the biosensor surface, the opposite of association and the average kdis of 380.5 μM, 38.05 μM, 3.805 μM, 0.3805 μM, and 0.03805 μM IVM were 8.20 × 10 -3 s −1 , 1.20 × 10 -5 s −1 , 2.38 × 10 -4 s −1 , 2.04 × 10 -4 s −1 , and 2.16 × 10 -4 s −1 , respectively.Figure 1c shows that the 380.5 μM IVM concentration demonstrates a significantly larger dissociation than all other concentrations (p < 0.0001).This difference is in line with what is seen in Fig. 1a.In the other studies mentioned above, the various concentrations of IVM ranging from 3.13 to 100 μM combined with their respective proteins demonstrated a binding curve that more closely resembles that of the 380.5 μM concentration 52 .Finally, the average equilibrium dissociation constants (KD) are determined by the system based on the ka and kdis values.For the IVM concentrations in decreasing concentration, the KD values were found to be 1.61 × 10 -4 M, 1.21 × 10 -12 M, Vol:.( 1234567890 ), the trend begins to slowly increase again for the decreasing IVM concentrations.KD is inversely related to the affinity, and the higher the affinity the stronger the attraction between the drug and receptor.There are various uses for drugs varying within the micromolar to picomolar range and larger affinities do not always mean a more successful drug.For example, a study from 1984 found that the micromolar over nanomolar affinity of benzodiazepine behaved as a Ca 2+ channel antagonist 57 .

Binding kinetics of RBD
The 1:1, 1:2, 1:5, and 1:10 dilutions of the original RBD were determined to have respective concentrations of 26.8 μM, 13.4 μM, 5.36 μM, and 2.68 μM based on the molecular weight of 58 kDa.The basic kinetics test was repeated with the RBD samples to determine a baseline on how RBD at these concentrations interacts with the hydrophobic sensors.At the beginning of the association step, all samples undergo significant binding before leveling out.By the end of the association step the average binding layers were found to be 2.170 nm, 2.256 nm, 2.115 nm, and 2.066 nm for 26.8 μM, 13.4 μM, 5.36 μM, and 2.68 μM RBD, respectively (Fig. 2a).This level of binding was overall much thicker than that of IVM seen in Fig. 1a, and stayed relatively unchanged over the testing period.The similarities in the binding curves are illustrated in Fig. 2b where the average ka values for 26.8 μM, 13.4 μM, 5.36 μM, and 2.68 μM RBD were found to be 3.10 × 10 5 Ms −1 , 4.41 × 10 5 Ms −1 , 2.80 × 10 5 Ms −1 , and 8.82 × 10 5 Ms −1 , respectively.The differences in association with varying RBD concentrations were not statistically significant (p = 0.124).On the other hand, a general decrease in dissociation with decreasing RBD concentration was found (p = 0.027).The corresponding kdis values for 26.8 μM, 13.4 μM, 5.36 μM, and 2.68 μM RBD were 9.65 × 10 -5 s −1 , 5.71 × 10 -5 s −1 , 2.45 × 10 -5 s −1 , and 1.00 × 10 -6 s −1 , respectively (Fig. 2c).Due to the similar values for the association values, the magnitude of the KD values was determined by the decreasing dissociation values with decreasing RBD concentrations (p = 0.023).Accordingly, average KD values were determined to be 3.78 × 10 -10 M, 1.54 × 10 -10 M, 9.33 × 10 -11 M, and 2.85 × 10 -12 M. Previous study performed on SARS-CoV RBD and hACE2 receptors used concentration ranges between 1.85 nM to 1.67 µM, overall smaller than the concentrations used in this experiment [58][59][60] .This difference is primarily seen in the thickness of the binding level which is much smaller than what is seen in this experiment.There is a consistent increase in binding demonstrated in all samples with increasing RBD concentration, the same general trend seen in this experiment.Additionally, the association rates are very similar with all results on the same magnitude of 10 5 Ms −1 .However, the dissociation values were consistently higher than the values seen in this experiment as they were consistently in the magnitude of 10 -3 s −1 .This is most likely due to the fact that the hACE2 was used instead of APS biosensors and other differences in the experimental setup.Due to the increased dissociation rates found in the studies, the KD values were also decreased to the nanomolar range.Overall, our findings are consistent with these studies with some differences due to the experimental setup.

Binding kinetics of ivermectin and RBD combined
To observe the effect of IVM on the binding behaviors between RBD and the APS sensors, the basic kinetics test was performed with a 50/50 combination of the two.RBD at 5.36 μM concentration was used to mix with each of the IVM concentrations.Comparing Figs.2a and 3a, it can be seen that overall the binding layers were slightly lower when the RBD was combined with IVM versus when it was RBD alone but still noticeably higher than IVM by itself.The average binding layers at the end of the association stage were 1.976 nm, 2.069 nm, 2.085 nm, 2.002 nm, and 2.039 nm for the 380.5 μM, 38.05 μM, 3.805 μM, 0.3805 μM, and 0.03805 μM IVM combined with 5.36 μM RBD, respectively.As seen in Fig. 3b, there is a consistent increase in association, with corresponding decreasing IVM concentrations (p = 0.0005).This trend is in line with what is seen in the IVMonly trials but demonstrates a more consistent increase.The ka values were found to be 4.12 × 10 3 Ms −1 , 6.90 × 10 4 Ms −1 , 4.74 × 10 5 Ms −1 , 5.03 × 10 6 Ms −1 , and 3.47 × 10 7 Ms −1 for the decreasing IVM concentrations combined with RBD.Compared to the RBD-only results the combination with the 3.805 μM IVM provides a minimal effect on the association between RBD and the APS biosensor.However, the smaller concentrations demonstrate an increase in association while the larger concentrations show a decrease.One way to reduce the efficacy of viral replication is through competitive inhibition that defers the association of the virus to the body 61,62 .This negative cooperativity is only seen at the higher concentrations 38.05 μM and 380.5 μM.For the concentrations lower Only at approximately 100 and 1,000 times the currently approved dosage are the desired effects found, which is in line with the findings of previous studies that the current dosage, and even 3 times higher, do not seem to have much of an effect 26,63 .In one study that tested IVM at 1200 μg/kg, 6 times the approved dosage, there was no significant reduction in viral load but more participants dropped out due to the tolerability of the high dosage and did not recommend additional trials at such a high level be performed 64 .While IVM is presumed to inhibit the nuclear translocation of viral proteins, clinical efficacy has not yet been proven 65,66 .The effect of IVM over an approximate 35-min pre-incubation time frame was determined using binding kinetics.Based on the pre-incubation time of each sample, the combination trial results were split into the following three groups: 5 min, 20 min, and 35 min.At the end of the association step the binding layers were between 1.842 nm and 2.175 nm for all concentrations and times (Supplementary Fig. 1).The statistical significance of the effect of IVM concentration or pre-incubation time on ka, kdis, and KD were obtained using ANOVA tests and listed in Supplementary Table 1.A general increasing trend of association for decreasing IVM concentrations was observed in the time trials (p = 0.0013-0.098).While, nominal changes in the association values of the IVM and RBD complex to the hydrophobic APS sensors over time for each given concentration (p = 0.144-0.67),as seen in Fig. 4.There was a correlation between decreasing IVM concentration and decreasing kdis values (p = 0.008-0.0858).The low kdis values for the 0.03805 μM IVM with 5-min pre-incubation could be explained by the hydrophobicity analysis from the MD simulation in the following sections.With a small amount of IVM, the hydrophobic residues that were initially embedded in the RBD were attracted outwards to interact with the IVM, which in turn increased the hydrophobic area on the RBD surface.The hydrophobic regions were then occupied by the IVM molecules when IVM concentration was sufficient, leading to a higher kdis value as shown with the higher IVM concentration results.When the IVM concentration was low and depleted, however, RBD binding to the hydrophobic sensor became stronger which resulted in a low kdis value.Additionally, no discernable trend was observed for the dissociation values of each IVM concentration combination with RBD over the time frame (p = 0.144-0.699).Despite the fluctuations of the association and dissociation values, KD remained relatively unchanged throughout the testing process (p = 0.291-0.798).Each of the average KD time trial results for the individual IVM concentration is within the same order of magnitude as that of the combined combination trials and following the same increasing trend with decreasing IVM concentration (p = 0.0036-0.015).Therefore, pre-incubation time did not have a statistically significant effect on ka, kdis, or KD at all tested IVM concentrations.One potential reason for the minimal changes found across the span of this experiment is the relatively short incubation period of 5-35 min compared to 46 h for the in vitro study and 0 to 24 h in animal studies 13,29,30 .However, in comparison to MD simulations of the interactions that are typically performed in the ns range, it is much longer 8,36 .It is important to consider the kinetics, interactions, and effects of IVM over a wide scope of time differences to ensure efficacy and safety when considering repurposing IVM for use with COVID-19.

Protein-ligand docking verification
The highest binding affinity of IVM with RBD was − 4.73 kcal/mol.The docking pose of the IVM to the RBD at the binding site was shown in 2D and 3D space.Multiple non-covalent bonds were formed between the IVM and RBD, including four hydrogen bonds at residues LYS 403, TYR 453, GLU 484, and SER 494, hydrophobic interactions at TYR 449, CYS 488, TYR 489, PHE 490, TYR 495, PHE 497, TYR 505, PHE 456, and LEU 455, and one polar bond at GLN 493 (Fig. 5a).The distances of the hydrogen bonds were 1.66 Å, 2.01 Å, 2.01 Å, and 1.87 Å at residues GLU 484, SER 494, LYS 403, and TYR 453, respectively (Fig. 5b).Our results agree with the study of Saha and Raihan 67 that LEU 455, a residue that favors interaction with human ACE2, can bind to IVM.The binding between RBD and IVM is largely regulated by hydrophobic interactions, similar to that of the human ACE2-IVM complex that was reported in a previous study 68 .Many of the residues-GLN 493, TYR 505, TYR 449, TYR 489, PHE 456, TYR 495, and LEU 455-that were found to regulate the binding between RBD and IVM were previously identified as hotspots in interactions between RBD and ACE2, which made IVM a candidate for drug repurposing for COVID-19 69 .

Molecular dynamics
All-atom molecular dynamics in an explicit solvent model with a total length of 100 ns was performed on three protein-ligand complexes to further understand the interaction between IVM and RBD.In order to determine which water model to be used in the full-length simulation, short MD simulations of 10 ns were performed on the docked RBD and IVM complex with three widely exploited explicit water models-SPC, SPC/E, and TIP3P.The protein Cα RMSD values with the three solvent models were compared and all three models showed similar RMSD values and fluctuation throughout the 10 ns simulations (Supplementary Fig. 2).Therefore, the TIP3P model was chosen for the full-scale 100 ns MD simulations as the water molecules in the TIP3P model have rigid geometry that closely matches the actual situation 70 .0.1 M explicit Na + and Cl − ions were added in the solvent model and two additional Cl − ions were added to neutralize the system.The MD simulations were performed on three complexes-unbound RBD, RBD with one IVM, and RBD with five IVM molecules-under a temperature of 300 K and pressure of 1 atm for 100 ns.The unbound RBD served as a control to compare for any structural changes of the protein that was induced by IVM.The initial and final snapshots of the superimposed structures of the RBD protein bound with IVM showed that IVM can bind and stay closely in the active site of the RBD protein (Fig. 6a,b).The binding energy between the RBD protein and the IVM was calculated using the MM-GBSA method from the last 40 ns of the MD simulation trajectory.The overall average binding free energy between the RBD protein and the IVM was − 50.20 ± 8.32 kcal/mol, which was considered as high binding energy by other studies and supported the protein-ligand docking result of high binding affinity between the RBD protein and the IVM [71][72][73] .The RBD protein Cα RMSD values of the three complexes were monitored throughout the MD simulations and analyzed to show the stability of protein structure (Fig. 6c).All simulations were performed in three replicates and the other two data sets of each system were included in Supplementary Fig. 3-5.The average RMSD value of the unbounded RBD protein, RBD protein with one IVM molecule, and RBD protein with five IVM molecules were 2.77 ± 0.44 Å, 2.48 ± 0.41 Å, and 2.74 ± 0.28 Å, respectively, and all the MD simulations were considered stable after 60 ns.In the unbound RBD, the RMSD value increased during the first 25 ns and reached a high value of 4 Å, then it reached equilibrium around 2.7 Å with small fluctuation until the end of the simulation.The complex of RBD and one IVM molecule maintained a stable RMSD of 2.2 Å for the first half of the simulation, then the RMSD value increased slightly after 55 ns up to around 3 Å.For RBD with five IVM molecules, the RMSD value increased during the first 20 ns and remained stable at around 3 Å for the rest of the simulation.Comparing the two complexes with IVM to the control, it was shown that the RMSD values of the three systems overlapped during the second 50 ns of the simulations.This indicates that the protein structure remained relatively the same during the simulations.However, RBD with five IVM took shorter time to stabilize than with one IVM molecule, suggesting that the increased IVM concentration in the system was beneficial for the protein stability.The protein Cα RMSF was calculated over the trajectory to investigate the structural fluctuation in terms of residues.The two protein-ligand complexes shared the same RMSF features as the unbound RBD alone, with residues 475 to 486 fluctuating the most during the simulations (Fig. 6d).The location and fraction of RBD interaction with the IVM molecule during the simulation were categorized into three interaction types-hydrogen bonds (H bonds), hydrophobic interaction, and water bridge (Fig. 6e).During the 100 ns simulation, the interaction between IVM and residue 493 was maintained the entire time and with residue 490 for 60% of the time.Other residues that maintained an interaction over 20% of the time were residues 405, 484, 489, 496, and 505.Residues 484 and 493 were shown to play critical roles in the RBD and human ACE2 binding 39 .The strong contact and stable interaction between these two residues and IVM indicated that IVM might be able to prevent and interrupt the RBD-ACE2 binding, hence that several published in silico studies reported IVM as a potential treatment for COVID-19 67,68,74,75 .However, more experimental evidence  www.nature.com/scientificreports/ is still required to further investigate the effect of IVM in preventing or treating SARS-CoV-2, and this study attempts to fill this gap.
To further study the structural properties of the RBD protein in the systems with or without IVM molecules, post-MD analysis such as SASA and Rg were performed.SASA of a protein is defined as the protein surface area that is accessible to the solvent and is a measurement of protein stability 76 .The average SASA values for the unbounded RBD protein, RBD protein with one IVM, and RBD protein with five IVM were 11,067.36± 210.72 Å 2 , 11,471.19 ± 216.84 Å 2 , and 13,085.14± 324.78 Å 2 , respectively.There was no major change of SASA values between the unbound RBD protein and the RBD protein with one IVM molecule, however, the SASA of the RBD protein with five IVM molecules was higher than the other two systems (Fig. 6f).The RBD protein with five IVM molecules having higher SASA therefore indicated a less stable overall protein structure compared to the unbound RBD protein and the RBD protein with only one IVM, while the RBD protein in the other two systems exhibited no major difference in protein stability.The Rg measures the compactness of the protein structure and is another indicator for protein stability.The average Rg values for the unbounded RBD protein, RBD protein with one IVM, and RBD protein with five IVM were 18.50 ± 0.14 Å, 18.74 ± 0.17 Å, and 19.52 ± 0.17 Å, respectively.Similar to the SASA, the RBD protein in the complex with five IVM molecules exhibited higher Rg compared to the other two systems, while no major difference in Rg was observed between the unbound RBD protein and the RBD protein with one IVM molecule.Although the RBD with five IVM molecules had an overall higher Rg, the Rg value showed low variation and stayed steady after 50 ns of the simulation, suggesting a still relatively stable protein structure (Fig. 6g).
Pairwise distance PCA was performed over the 100 ns trajectory to further study the motion of the protein system in the MD simulations (Fig. 7a-c).PCA has been used to describe the total fluctuations of systems during MD simulation with drastic dimensionality reduction [76][77][78] .The first two eigenvectors PC1 and PC2 accounted for ~ 40% of the total variance in the three systems-unbound RBD protein, RBD protein with one IVM molecule, and RBD protein with five IVM molecules-and were used to construct the PCA.The PCA plots of all the three systems demonstrated smooth and consecutive changes of the data points with no outliers, with the clustering of data points occurring after 60 ns (light green and yellow points) suggesting stabilization of the protein structure.The concentrated cluster distribution in the last 40 ns of the RBD protein with one IVM was an indicator of reduced protein flexibility possibly due to ligand binding, which were in agreement with the other parameters such as RMSD and SASA.
The initial and final frame of the MD simulation of RBD with five IVM molecules were analyzed for the protein surface hydrophobicity.The IVM molecules moved towards RBD and attached to the hydrophobic RBD surface at the end of the simulation (Fig. 8a,b).The area of the hydrophobic RBD surface increased from 846.5 Å 2 at the initial frame to 1143.5 Å 2 at the final frame when interacting with IVM molecules.It was hypothesized that the increased hydrophobicity of RBD was due to IVM attracting the embedded hydrophobic residues to come out to interact with the ligand.By including one and more than one (five) IVM molecules in the system, we aimed to study the effect of the increasing IVM concentration on the protein stability and surface properties which might explain the different kinetic parameters that were experimentally observed.A maximum of five IVM molecules were tested in this study due to the constraints of time and resources.With more IVM in the solvent system, it is likely that a majority of the RBD hydrophobic area would be occupied by the ligands, which in turn reduces the binding ability of the RBD to the APS biosensor.This could explain why increasing IVM concentration led to decreasing rate of association of RBD in the experiment.It is important to note that, with increasing numbers of IVM molecules in the in silico system, the equivalent dosage in human body will exceed the established dosage limit of 150-200 μg/kg that was proved to have adverse effects for human 33 .

Conclusions
The binding behaviors are an instrumental aspect when developing effective treatment or prevention methods for the virus.This study utilized BLI technology to determine the binding kinetics and affinity of various IVM concentrations to a hydrophobic biosensor and their effect when combined with RBD.Negative cooperativity, where IVM slows or prevents RBD from binding to the biosensor, was only found at 38.05 μM and especially at 380.5 μM, approximately 100 and 1,000 times the approved IVM dosage.These same concentrations were able to weaken the affinity levels between RBD and the hydrophobic biosensor as desired, while the lower concentrations had a minimal effect or strengthened it.Only nominal effects on the attachment and detachment of the IVM and RBD complex to and from the hydrophobic APS sensors were seen due to changes in pre-incubation time.The strong binding capability between IVM and RBD was verified through in silico molecular docking.The RBD residues that interacted with IVM were found to be hotspot residues in the RBD-ACE2 complex, showing that IVM might be able to interrupt the SARS-CoV-2 infection in human cells which makes IVM a potential prevention and treatment of COVID-19.MD simulation of the unbound RBD and RBD with one and five IVM molecules demonstrated good stability of the complex structure over the 100 ns simulation time.Comparison analysis between the three complexes further revealed that higher IVM concentration triggers an increase in the hydrophobic surface area of RBD, which changes its binding ability.Overall, based on our results, it is not recommended that IVM, at the currently approved dosage, be used as a preventative or treatment method.The established dosage used in this experiment is based on oral administration.This provides the drug with the opportunity to disperse throughout the body, while COVID-19 on the other hand is centralized in the respiratory system.Other methods, such as nasal spray or nebulizer treatment, should be considered as they are better suited for direct delivery to the respiratory system.Potential future studies could investigate IVM on the different concentrations of RBD to mimic the varying viral loads.This change, in addition to more extended time trials, will allow a further understanding of the effects of IVM throughout the various stages of the viral infection.
https://doi.org/10.1038/s41598-024-53086-0www.nature.com/scientificreports/7.51 × 10 -11 M, 1.86 × 10 -10 M, and 4.81 × 10 -11 M, respectively.The combination of the lowered association and increased dissociation of the 380.5 μM concentration leads to a much higher KD value, in the micromolar range, compared to the other concentrations in the picomolar range.While there is a large decrease in the KD values between the 380.5 and 38.05 μM concentrations (p = 0.0003

Figure 1 .
Figure 1.The binding kinetics for IVM only attach to and detach from the hydrophobic APS sensors.(a) The binding curves of IVM were generated using the basic kinetic analysis of IVM only at 380.5 μM, 38.05 μM, 3.805 μM, 0.3805 μM, 0.03805 μM.The vertical dotted line indicates the transition from the association step to the dissociation step.In the same order of the three groups, the association constant ka is shown in (b).The dissociation constant kdis is shown similarly in (c).

Figure 2 .
Figure 2. The binding kinetics for RBD only attach to and detach from the hydrophobic APS sensors.(a) The binding curves of IVM were generated using the basic kinetic analysis of RBD only at 26.8 μM, 13.4 μM, 5.36 μM, and 2.68 μM.The vertical dotted line indicates the transition from the association step to the dissociation step.In the same order as the three groups, the association constant ka is shown in (b).The dissociation constant kdis is shown similarly in (c).

Figure 3 .
Figure 3.The binding kinetics for IVM and RBD combination attachment to and from the hydrophobic APS sensors.(a) The binding curves of IVM were generated using the basic kinetic analysis of 80.5 μM, 38.05 μM, 3.805 μM, 0.3805 μM, and 0.03805 μM IVM combined with 5.36 μM RBD.The vertical dotted line indicates the transition from the association step to the dissociation step.In the same order as the three groups, the association constant ka is shown in (b).The dissociation constant kdis is shown similarly in (c).

Figure 4 .
Figure 4.The effects of pre-incubation time on the attachment and detachment of the IVM and RBD complex to and from the hydrophobic APS sensors.(a) The association constant ka and (b) the dissociation constant kdis of 5.36 μM RBD and 380.5 μM, 38.05 μM, 3.805 μM, 0.3805 μM, and 0.03805 μM IVM combined with 5 min, 20 min, and 35 min pre-incubation.

Figure 5 .
Figure 5.The docking pose of IVM with SARS-CoV-2 RBD in (a) 2D and (b) 3D interaction diagrams.The hydrogen bonds formed between IVM and RBD residues from the left to right in the 3D interaction diagram are SER 494, LYS 403, TYR 453, and GLU 484, and the corresponding bond lengths are shown in purple.The RBD residues that interacted with the IVM are labeled.

Figure 6 .
Figure 6.Structural properties of the protein-ligand complexes during MD simulation.(a) Initial and (b) final snapshots of RBD protein bound with IVM.(c) The RBD protein Cα RMSD values over time for the average displacement change in the unbound RBD (red), RBD with one IVM molecule (yellow), and RBD with five IVM molecules (blue).(d) Residue-based protein Cα RMSF over the trajectory for local fluctuation along the protein chain in the unbound RBD (red), RBD with one IVM molecule (yellow), and RBD with five IVM molecules (blue).(e) Fraction of RBD interactions with the IVM molecule over the trajectory categorized into hydrogen bonds (green), hydrophobic interactions (orange), and water bridges (blue).(f) SASA and (g) Rg of the RBD protein over time in the unbound RBD (red), RBD with one IVM molecule (yellow), and RBD with five IVM molecules (blue).

Figure 7 .
Figure 7. Pairwise distance PCA on the RBD protein Cα from the MD simulation trajectory.The first two eigenvectors PC1 and PC2 were used to construct the PCA in (a) unbound RBD, (b) RBD and one IVM molecule, and (c) RBD and five IVM molecules.

Figure 8 .
Figure 8. Hydrophobic (orange) and hydrophilic (cyan) surface of RBD interacting with 5 IVM molecules (green sticks) in the (a) initial and (b) final frame of MD simulation.